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ABSTRACT 

This  work  is  concerned  with  hyperbolic  systems  of  partial  differential 
equations  for  which  certain  of  the  associated  propagation  speeds  are  a  great 
deal  larger  than  the  other  propagation  speeds.  Our  goal  is  to  find  boundary 
conditions  which  prevent  rapidly  moving  waves  from  entering  the  given  spatial 
domain.  Conditions  of  this  type  are  desirable  in  certain  numerical 
computations  arising  in  meteorology. 

In  order  to  find  these  conditions,  we  first  transform  the  given  system  to 
an  approximate  diagonal  form  in  such  a  way  that  each  of  the  new  dependent 
variables  can  be  identified  as  a  slow,  incoming  fast,  or  outgoing  fast  com¬ 
ponent  of  the  solution.  We  then  find  local  boundary  conditions  which  suppress 
the  incoming  fast  part.  Pseudo-differential  operators  are  used  throughout  the 
entire  process.  We  consider  only  linear  systems.  In  Part  II  of  this  work  [7] 
these  methods  are  applied  in  detail  to  the  linearized  shallow  water  equations. 

The  results  of  numerical  tests  of  various  boundary  conditions  are 
included  in  both  papers.  We  also  outline  a  method  by  which  the  conditions  can 
be  justified  analytically. 
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SIGNIFICANCE  AND  EXPLANATION 


When  one  computes  a  numerical  approximation  to  the  solution  of  a  partial 
differential  equation,  it  is  sometimes  .acessary  to  restrict  the  computation 
to  a  domain  which  is  only  a  portion  of  the  domain  on  which  the  problem  is 
naturally  defined.  This  is  done  when  one  is  interested  in  the  behavior  of  the 
solution  only  on  that  part  and  when  a  computation  over  the  entire  domain  would 
be  prohibitively  expensive. 

In  such  cases  a  portion  of  the  boundary  of  the  computational  domain 
represents  merely  the  edge  of  the  computation  and  corresponds  to  nothing 
physical.  The  task  of  finding  suitable  boundary  conditions  to  impose  on  this 
part  of  the  boundary  can  present  some  analytical  and  numerical  difficulties. 
This  work  is  concerned  with  one  such  difficulty  which  arises  in  limited-area 
meteorological  computations.  It  is  also  related  to  the  problem  of  finding 
conditions  which  prevent  the  nonphysical  reflection  of  waves  at  an  artificial 


boundary. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  author  of  this  report. 


BOUNDARY  CONDITIONS  FOR  SUPPRESSING  RAPIDLY 
MOVING  COMPONENTS  IN  HYPERBOLIC  SYSTEMS,  I 


Robert  L.  Higdon 


1.  INTRODUCTION 

Hyperbolic  partial  differential  equations  are  characterized  by  the  fact 
that  they  propagate  information  at  finite  speed.  For  first  order  hyperbolic 
systems  there  may  be  several  such  propagation  speeds,  each  corresponding  to  an 
eigenvalue  of  the  principal  symbol  of  the  system.  In  the  present  work  we  con¬ 
sider  systems  for  which  the  various  speeds  can  have  substantially  different 
magnitudes.  Systems  of  this  type  are  sometimes  said  to  have  "multiple  time 
scales". 

Examples  of  such  systems  arise  in  the  study  of  fluid  dynamics.  The 
shallow  water  equations  and  the  Euler  equations  of  gas  dynamics  admit  certain 
modes  associated  with  the  movement  of  the  fluid  and  certain  other  modes 
associated  with  the  propagation  of  gravity  waves  and  sound  waves,  res¬ 
pectively.  If  these  waves  move  at  speeds  which  are  considerably  greater  than 
rate  of  flow  of  the  fluid,  then  these  systems  have  two  time  scales. 

In  this  paper  we  consider  initial-boundkry  value  problems  for  such 
systems.  The  goal  is  to  find  boundary  conditions  which  prevent  rapidly  moving 
high-frequency  waves  from  entering  the  given  spatial  domain.  That  is,  we  will 
try  to  identify  the  portion  of  the  solution  which  is  entering  the  region  at 
the  fast  speed,  and  we  will  then  attempt  to  set  this  part  of  the  solution 
equal  to  zero  at  the  boundary.  The  purpose  of  these  conditions  is  to  prevent 
excessive  errors  in  the  boundary  data  from  propagating  rapidly  into  the 
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interior  during  numerical  computations*  In  Section  2  of  this  paper  we  des¬ 
cribe  a  physical  problem  in  which  boundary  conditions  of  this  type  would  be 
useful/  and  at  the  end  of  that  section  we  mention  a  connection  between  this 
work  and  the  construction  of  nonreflecting  boundary  conditions  for  hyperbolic 
systems . 

In  order  to  find  the  conditions  desired  here/  we  first  transform  the 
given  system  to  an  approximate  diagonal  form  so  that  at  high  frequencies  each 
of  the  new  dependent  variables  can  be  identified  as  a  slow,  incoming  fast,  or 
outgoing  fast  component  of  the  solution.  Certain  portions  of  this  uncoupling 
process  are  similar  to  some  ideas  outlined  by  Engquist  and  Majda  [1]  for 
constructing  absorbing  boundary  conditions  for  first  order  systems.  After  the 
system  is  uncoupled  we  find  local  boundary  conditions  which  suppress  the 
incoming  fast  part.  Pseudo-differential  operators  are  used  extensively 
throughout  the  entire  process.  The  discussion  is  limited  to  linear  systems. 

In  Sections  3  and  5  of  this  paper  we  consider  problems  in  one  and  several 
space  dimensions,  respectively.  In  order  to  simplify  the  discussions  in  these 
sections,  we  avoid  using  pseudo-differential  operators  and  instead  simulate 
their  use  with  formal  manipulations yof  Fourier  transforms.  The  calculations 
with  Fourier  transforms  give  results  which  are  valid  for  systems  having  con¬ 
stant  coefficients,  aside  from  a  technical  difficulty  about  the  use  of  the 
transforms  which  is  mentioned  in  Section  3.2.  Pseudo-differential  operators 
are  used  to  settle  this  point  and  to  generalize  the  methods  of  Sections  3  and 
5  to  systems  which  have  variable  coefficients.  The  generalization  is  outlined 
briefly  in  Section  4. 1  and  is  developed  in  detail  in  Part  II  of  the  this  work 
[7].  There  we  assemble  the  ideas  of  Section  5  and  give  a  detailed  development 
of  boundary  conditions  for  the  linearized  shallow  water  equations,  or, 
equivalently,  the  two-dimensional  Euler  equations  of  gas  dynamics. 


The  effects  of  various  boundary  conditions  derived  here  can  be  analyzed 
using  techniques  from  the  theory  of  propagation  of  singularities  for  linear 
partial  differential  equations.  This  analysis  is  outlined  in  Section  4.2  of 
this  paper,  and  it  is  discussed  in  greater  detail  in  the  thesis  [6]  on  which 
these  two  papers  are  based.  We  also  present  the  results  of  some  numerical 
tests  cf  the  boundary  conditions,  in  Section  3.4  for  a  model  problem  in  one 
space  dimension  and  in  Part  XI  for  the  shallow  water  equations. 


This  work  is  motivated  primarily  by  a  situation  which  can  arise  when  one 
computes  numerical  approximations  to  the  solutions  of  hyperbolic  systems  which 
model  the  behavior  of  the  earth's  atmosphere.  Under  certain  circumstances  the 
boundary  data  available  for  these  computations  are  substantially  less  accurate 
than  the  initial  data  which  are  available.  The  reason  for  this  will  be 
discussed  below.  This  situation  is  unfortunate,  since  the  large  errors  in  the 
boundary  data  can  generate  comparable  errors  in  the  interior  and  thereby  waste 
the  extra  accuracy  contained  in  the  initial  data.  The  goal  of  the  present 
work  is  to  reduce  as  much  as  possible  the  extent  to  which  this  contamination 
takes  place. 

We  propose  to  do  this  by  imposing  boundary  conditions  which  prevent 
rapidly  moving  waves  from  propagating  into  the  interior.  The  fast  waves 
allowed  by  the  systems  in  question  are  sound  waves  or  gravity  waves,  and  in 
the  physical  problem  these  waves  represent  only  a  small  portion  of  the  exact 
solution.  But  the  fast  modes  are  still  able  to  carry  substantial  errors  into 
the  spatial  domain,  since  the  errors  in  the  boundary  data  are  not  related  to 
the  exact  solution  and  in  general  can  excite  all  of  the  modes  in  the  system. 

We  would  therefore  like  to  identify  the  portion  of  the  solution  which  is 
entering  the  domain  at  the  fast  speed  and  then  set  this  portion  equal  to  zero 
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at  the  boundary.  This  would  have  the  desired  numerical  effects  and  would  also 
be  physically  realistic. 

We  will  not  attempt  to  prevent  the  propagation  of  error  at  the  slow 
speed,  since  it  is  not  realistic  to  assume  that  the  slow  part  of  the  solution 
is  equal  to  zero.  However,  it  is  still  worthwhile  to  stop  the  rapidly  moving 
portion,  since  in  meteorological  problems  the  fast  and  slow  speeds  can  easily 
differ  by  a  factor  of  five  or  ten  to  one. 


The  large  errors  in  the  boundary  data  are  commonly  found  in  limited  area 
computations  which  are  used  to  predict  local  atmospheric  phenomena.  In 
meteorological  calculations  a  natural  domain  of  definition  for  initial- 
boundary  value  problems  is  the  entire  earth's  atmosphere.  However,  the  size 
of  this  domain  and  the  limitations  of  present-day  computing  machines  require 
that  global  computations  be  done  on  meshes  which  are  quite  coarse.  In  current 
practice  the  grid  spacing  for  such  calculations  is  often  one  interval  per  two 
and  a  half  degrees  latitude  and  longitude.  These  calculations  can  give  useful 
information  about  global  phenomena,  but  the  grid  spacing  is  too  coarse  for 
predicting  local  phenomena. 

It  is  common  practice  to  perform  additional  computations  over  smaller 
regions  with  finer  meshes  so  that  these  local  phenomena  can  be  resolved.  For 
such  a  computation  the  spatial  region  is  a  cylinder  in  the  atmosphere  which  is 
bounded  by  the  earth's  surface,  the  top  of  the  atmosphere,  and  an  artificial 
computational  boundary.  The  artificial  boundary  merely  defines  the  edge  of 
the  computation  and  represents  nothing  physical.  The  problem  considered  here 
arises  when  one  tries  to  find  suitable  conditions  to  impose  on  this  additional 
portion  of  the  boundary.  On  this  section  it  is  necessary  to  prescribe  values 
for  variables  which  in  some  sense  represent  flow  into  the  region,  and  for  this 
one  needs  values  of  the  solution  on  the  boundary  for  times  occurring  after  the 
initial  time.  This  information  is  interpolated  from  the  results  of  the  global 
computation.  Unfortunately,  these  results  contain  large  errors  resulting  from 
the  use  of  the  coarse  mesh,  and  they  may  also  be  affected  by  inadequate 
initial  data  on  portions  of  the  atmosphere.  On  the  other  hand,  the  initial 
data  on  the  computational  domain  can  be  considered  fairly  reliable,  since  one 
would  perform  local  computations  only  over  a  populated  region  where  there  is  a 
dense  network  of  observation  stations  capable  of  making  accurate  measurements. 


This  explains  the  substantial  difference  between  the  quality  of  the  initial 
data  and  the  quality  of  the  boundary  data. 

The  fast  modes  allowed  by  the  system  can  cause  another  numerical  dif¬ 
ficulty  besides  the  one  discussed  above.  Because  of  the  Courant-Friedrichs- 
Lewy  condition,  the  fast  modes  can  impose  a  severe  restriction  on  the  per¬ 
missible  time  step  for  stable  explicit  difference  approximations  to  the 
differential  equation.  In  general,  this  presents  the  choice  of  either  using 
an  implicit  difference  method  or  an  explicit  method  with  very  short  time 
steps*  Both  choices  are  undesirable  because  of  the  computational  expense 
which  is  involved.  However,  the  physical  insignificance  of  the  fast  waves 
makes  it  possible  to  use  a  certain  reduced  system  which  can  be  solved  much 
more  efficiently  than  the  original  system.  A  portion  of  the  reduced  system  is 
elliptic  and  must  be  solved  globally,  but  the  remainder  is  hyperbolic  with 
only  slow  propagation  speeds.  The  use  of  such  a  system  is  based  on  a  theory 
developed  by  Kreiss  [9]  for  problems  having  multiple  time  scales,  and  it 
requires  that  the  fast  scale  be  properly  filtered  from  the  initial  data.  The 
implementation  of  this  process  is  discussed  by  Gustafsson  and  Kreiss  [5]  and 
Gustafsson  [3] ,  [4] . 

This  method  does  not  quite  solve  the  problem  discussed  in  this  paper, 
since  the  elliptic  part  of  the  reduced  system  can  cause  certain  portions  of 
the  errors  in  the  boundary  data  to  propagate  instantaneously  into  the 
interior.  In  this  paper  we  will  be  concerned  mainly  with  the  high-frequency 
portion  of  the  solution,  since  numerical  errors  are  mainly  of  this  type. 

Errors  having  high  frequencies  in  the  tangential  variables  will  be  damped  by 
the  elliptic  part  of  the  system,  but  this  is  not  the  case  for  errors  in  t. 

The  methods  of  this  paper  are  oriented  mainly  toward  high  frequencies  in  time. 


It  might  seem  that  we  could  best  deal  with  the  problems  described  in  this 
section  by  modifying  the  system  of  differential  equations  so  as  to  prohibit 
solutions  containing  rapidly  moving  waves.  However,  it  is  not  known  how  to  do 
this  and  still  retain  a  system  which  is  mathematically  well-behaved  and  which 
is  a  sufficiently  accurate  model  of  the  atmosphere  to  be  useful  for 
meteorological  computations . 

He  conclude  by  noting  that  the  work  in  this  paper  is  relevant  to  the 
problem  of  finding  boundary  conditions  which  prevent  the  nonphysical 
reflection  of  waves  by  an  artificial  boundary.  A  wave  leaving  the  spatial 
domain  will  be  reflected  back  into  the  interior  if  and  only  if  the  boundary 
conditions  allow  the  outgoing  modes  to  stimulate  the  incoming  modes.  The 
linkage  between  the  modes  is  broken  if  the  system  is  uncoupled  and  the 
incoming  modes  are  set  equal  to  zero  at  the  boundary.  Conditions  of  this  type 
are  realistic  if  the  solution  consists  entirely  of  outgoing  waves. 


SYSTEMS  IN  ONE  SPACE  DIMENSION 


3. 1  General  Remarks 

We  will  consider  the  hyperbolic  system 


(3.1) 


for  0<x<1,t>0.  This  can  also  be  written  wfc  =  Awx  +  Cw,  where 
w  =  (u,v)  e  R  .  The  entries  in  A  and  C  are  functions  of  x  and  t. 

In  order  to  simplify  the  notation  we  have  chosen  a  system  having  two 
scalar  components.  It  can  be  seen  easily  that  the  ideas  presented  in  this 
section  work  equally  well  for  systems  having  several  components. 

There  is  no  loss  of  generality  in  assuming  that  A  is  diagonal.  The 
system  is  hyperbolic,  so  A  has  real  eigenvalues  and  a  complete  set  of 
eigenvectors.  If  A  is  not  diagonal,  then  a  suitable  similarity 
transformation  and  change  of  dependent  variables  can  be  made  to  bring  the 
system  to  diagonal  form. 

We  assume  | a | < < | b |  and  a  <  0,  h  <  0.  The  first  assumption 
guarantees  the  presence  of  propagation  speeds  having  substantially 
different  magnitudes.  The  second  assumption  is  made  for  the  sake  of 
definiteness.  It  also  contains  the  assumption  that  det  A  j*  0,  i.e.,  that 
the  boundary  is  noncharacteristic. 

The  problem  is  to  identify  the  "fast"  part  of  the  solution  of  the 
system  and  then  find  boundary  conditions  which  suppress  this  as  much  as 
possible.  To  some  degree  this  can  be  done  via  the  method  of  charac¬ 
teristics.  If  C  were  diagonal,  then  the  system  (3.1)  would  consist  of 


two  independent  equations  for  u  and  v.  These  are  ordinary  differential 
equations  along  the  characteristics  curves  ^  ■  -a  and  ~  =  -b, 
respectively.  Initial  values  for  these  equations  are  provided  by  initial 
data  (t  -  0)  and  boundary  data  (x  *  0,  in  this  case)  for  (3.1). 

Since  | b | >> | a | ,  v  is  the  component  which  propagates  information  at  the 
fast  speed,  and  it  would  therefore  suffice  to  set  v  ■  0  at  the  boundary. 

However,  if  C  is  not  diagonal,  then  this  boundary  condition  may  not 
be  adequate.  In  this  case  u  can  act  as  a  forcing  function  in  the 
equation  for  v,  and  in  general  u  is  nonzero.  The  boundary  data  could 
thus  influence  the  solution  in  the  interior  at  the  fast  speed  by  first 
influencing  u,  which  in  turn  would  force  v. 

We  would  therefore  like  to  have  boundary  conditions  which  are  more 
effective  than  the  one  mentioned  above.  These  conditions  can  be  obtained 
by  transforming  the  system  so  as  reduce  the  coupling  found  in  the  lower 
order  term.  This  can  be  thought  of  as  a  process  of  identifying  more 
precisely  the  quantity  which  moves  slowly  and  the  quantity  which  moves 
rapidly.  Refined  boundary  conditions  can  then  be  derived  by  setting  the 
new  fast  variable  equal  to  zero  whenever  permissible  and  then  expressing 
this  condition  in  terms  of  the  original  unknowns  u  and  v.  In  the  next 
section  we  discuss  a  method  for  reducing  the  coupling  in  the  system,  and 
in  Section  3.3  we  produce  boundary  conditions  from  the  results  of  the 
uncoupling. 

3. 2  A  Method  for  Reducing  the  Oaupling  in  Lower-Order  Terms 

This  method  was  used  by  Taylor  [12]  to  reduce  the  coupling  in  systems 
of  pseudo-differential  equations  in  which  the  leading  symbol  is  already 
known  to  be  in  block-diagonal  form.  It  is  essentially  a  simple 


perturbation  argument,  valid  £or  large  frequencies,  which  is  disguised  by 
the  language  of  pseudo-differential  operators.  It  is  similar  to  some 
uncoupling  methods  used  by  Kreiss  [8]  and  O'Malley  and  Anderson  [11].  In 
this  sub-section  we  outline  the  technique  using  Fourier  transforms  in  a 
formal  manner.  A  general  version  of  the  perturbation  method  is  given  in 
the  Appendix. 

Before  this  method  can  be  used,  the  leading  order  part  of  the  system 
must  be  brought  to  diagonal  form,  or  at  least  to  a  suitable  block-diagonal 
form.  For  systems  in  one  space  dimension  this  step  is  trivial,  and  it  has 
already  been  accomplished  for  the  system  (3.1).  In  several  dimensions 
this  step  is  rather  complicated.  It  will  be  discussed  at  length  in 
Section  5. 

The  system  considered  here  is  wt  ”  Awx  +  Cw.  He  will  first  apply  a 
Fourier  transform  in  order  to  reduce  it  to  a  system  of  ordinary  dif¬ 
ferential  equations,  and  we  will  then  apply  the  perturbation  argument.  In 
this  section  we  assume  that  the  system  has  constant  coefficients. 

We  will  not  apply  a  transform  in  x  because  this  would  require 
information  about  the  solution  away  from  the  boundary.  That  would  not  be 
appropriate  in  a  discussion  of  boundary  conditions.  Instead,  we  will  use 
a  Fourier  transform  in  t.  The  use  of  such  a  transform  is  complicated  by 
the  fact  that  the  solution  can  grow  too  rapidly  as  t  ♦  *  for  the 
transform  to  be  defined  in  the  usual  sense,  but  it  is  possible  to  remedy 
this  problem  with  a  localization  argument  which  uses  properties  of 
pseudo-differential  operators.  Details  are  given  in  [6].  The  growth 
problem  could  also  be  solved  by  using  a  Laplace  transform,  but  we  will  not 
do  this  because  the  time  derivative  would  cause  the  initial  data  to  be 
introduced  into  the  transformed  equations.  That  would  be  an  undesirable 


complication.  It  is  convenient  to  use  pseudo-differential  operators  to 
handle  this  problem  since  we  were  going  to  introduce  these  operators, 
anyway,  in  order  to  treat  systems  having  variable  coefficients. 

The  system  (3.1)  can  be  written  as 

(3.2)  wx  =  A"1wt  -  A- 1 Cw 


In  terms  of  components  this  is 


(3.3) 


Introduce  formally  the  Fourier  transform  in  t.  Let  €  be  the  dual 
variable,  and  let  u,v,w  be  the  transforms  of  u,v,w.  If  the 
coefficients  in  the  equation  are  constant,  (3.2)  becomes 


(3.4) 


wx(x,5)  *  i?A  'w  -  A 

-  iC (A-1  -  j£A_1C)w 

=  i£R(i€)w. 


When  £  is  large  the  matrix  R(i£)  is  a  perturbation  of  the 
diagonal  matrix  A~^.  We  will  find  a  similarity  transformation  which 
reduces  the  coupling  caused  by  the  off-diagonal  elements,  and  we  will  then 
use  this  transformation  to  make  a  change  of  dependent  variable. 

Let  Q(i?)  =  I  +  (i€)”^M,  where  M  is  a  matrix  to  be  determined. 

For  large  ?,  Q  1  exists  and  has  the  expansion 

Q-1  *  I  -  +  o(-^). 

£ 
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Using  (3.4 ) ,  we  have 


(3.5) 


QRQ-1  -  (I  +  -  ~A-1C)(I  -  -gH  +  OiC2)) 

-  A*1  +  -j(MA“1  -  A_1M  -  A-1C)  +  0(5~2) . 


The  coupling  is  of  order  -2  if  MA" 
i.e.,  if 


A~^M  -  A-*C  is  diagonal. 


11 

m1 2 

21 

m22i 

*11  “12 
*21  ®22 


a  cn  a  c12 
b  1  c 21  b  'c22 


is  a  diagonal  matrix.  He  therefore  set  to  zero  the  off-diagonal  entries 
in  (3.6). 

(3.7)  m12b_1  -  a_1m|2  -  a-1c12  “  0  (row  1,  column  2) 

m21a_1  -  b_1m21  -  b^c-^  ■  0  ( row  2,  column  1) 

This  yields 


m  ■  ■  —  — ■ 

12  -1  -1 
b  -  a 


21  -1  -1 
a  -  b 


No  conditions  are  imposed  on  the  values  of  m^  and  m22*  For 
convenience  we  take  these  to  be  zero.  The  matrix  Q  is  then  given  by 


Q  -  I  +  ( iC )  M,  or 


51 


Equation  (3.5)  becomes 


(3.9) 


QRQ_1  -  q(a_1  -  (i5)“1A"1c)Q"1 


+ 


_i  /'r'cn 
15  V  0 


We  now  use  the  similarity  transformation  to  reduce  the  coupling  in 
equation  (3.4).  This  system  can  be  written  in  the  form 

(3.10)  (qw)  -  Q( 15a” 1  -  a“1c)q“1(qw) 

If  we  let  w^>  Qw  and  use  (3.9),  the  system  (3.10)  becomes 


When  €  is  large,  the  coupling  caused  by  the  lower  order  term  in 
(3.11)  is  weaker  than  the  coupling  in  the  original  system  (3.4).  This 
means  that  for  large  5  we  can  identify  more  precisely  the  rapidly  moving 
part  of  the  solution  and  do  a  better  job  of  suppressing  it.  This  will  be 
done  in  the  next  sub-section.  The  restriction  to  high  frequencies  is  not 
a  serious  one,  since  the  goal  of  this  work  is  to  suppress  the  effect  of 
numerical  errors,  which  are  mainly  high  frequency  phenomena. 

This  method  can  be  applied  repeatedly  to  reduce  even  further  the 
coupling  at  high  frequencies.  To  reduce  the  coupling  to  o(5  2) ,  we 
would  multiply  (3.11)  by  a  matrix  of  the  form  I  +  (i£)~Sl  and  then 
determine  M2  in  the  same  way  that  we  found  the  matrix  M  above.  In 
general,  to  reduce  the  coupling  from  o(5~n+1)  to  o(5  °),  we  would  use 


a  multiplier  of  the  form  I  +  5  nMn  The  details  of  this  process  involve 
no  new  ideas  and  will  not  be  given  here. 

The  method  has  been  presented  for  2x2  matrices.  It  can  also  be 
used  to  reduce  the  coupling  in  square  matrices  of  any  order  which  are 
perturbations  of  block-diagonal  matrices.  In  this  case  the  equations 
analogous  to  (3.7)  can  be  solved  provided  the  diagonal  blocks 
corresponding  to  a~^  and  b~^  have  disjoint  spectra.  This 
generalization  will  be  discussed  in  the  Appendix. 


3.3  Boundary  Conditions 

We  now  use  the  results  of  the  uncoupling  to  find  boundary  conditions 
which  suppress  the  fast  part  of  the  solution.  The  new  dependent  variable 
is  defined  by  w^  ■  Qw.  By  (3.8)  this  is 


(3.12) 


/u,(xf£)\ 

<«,«)/ 


/ 


\ 


\ 


_L  r! fji_ > 

U  '“b  -  a' 


/ 


I  u(x,£)1 
yvfx,^ 


v1  is  our  new  notion  of  what  constitutes  the  rapidly  moving  part  of  the 
solution.  For  large  €  it  is  a  perturbation  of  the  fast  characteristic 
variable  v.  To  suppress  the  fest  part  we  set  v^  ■  0  at  x  -  0,  i.e.. 


(3.13) 


AC 

1  r  21  ^  „ 

_  (— j)u  ♦  v  -  0  at 


To  obtain  local  boundary  conditions  we  multiply  (3.13)  by  iC  and  then 


apply  an  inverse  Fourier  transform.  The  result  is 


(3.14) 


It  +  “ 0  at  x  ‘ 0 


When  the  coefficients  in  the  system  depend  on  t,  this  derivation  is 
obviously  not  valid.  However,  it  is  possible  to  use  pseudo-differential 
operators  to  show  that  the  boundary  condition  (3.14)  still  has  some 
desirable  effects.  This  is  done  in  [6],  and  in  Part  II  [7]  we  do 
something  analogous  for  the  linearized  shallow  water  equations. 

In  order  to  define  a  well-posed  problem,  we  must  also  prescribe  a 
boundary  condition  for  the  slow  part  of  the  solution  of  (3.1).  In  this 
case  it  is  necesary  to  prescribe  a  value  at  x  *  0.  One  possible 
condition  is 


u  “  given  function  ,  for  x  ■  0. 

We  might  also  use  the  result  (3.12)  of  the  uncoupling.  That  is,  we  could 
prescribe  a  value  for  u^,  clear  denominators,  and  then  obtain  a  boundary 
condition  analogous  to  (3.14).  However,  there  seems  to  be  little  value  in 
doing  this,  so  this  possibility  will  be  disregarded. 

We  therefore  propose  the  boundary  conditions 


ac. 


<3‘15)  3?+(b-^u  "  ° 

u  m  g,  for  x  *  0, 

where  g  is  a  given  function  of  t.  An  integration  of  the  first  equation 
shows  that  v(0,t)  is  defined  in  terms  of  g  and  the  initial  value 
v(0,0).  The  conditions  (3.15)  thus  prescribe  values  for  the 
characteristic  variables  u,v  on  the  portion  of  the  boundary  where  the 
characteristics  enter  the  region.  This  implies  that  the  initial-boundary 
value  problem  is  well-posed. 


-15- 


3.4  Numerical  Computations 


In  this  sub-section  we  present  the  results  of  some  numerical 
computations  which  indicate  the  value  of  reducing  the  coupling  caused  by 
lower  order  terms.  He  consider  the  system 

«»•«>  nO  -  c-1  -.hSq?  ♦  SC) 

for  0  <  x  <  1  and  t  >  0.  This  is  the  system  (3.1),  where  a  ■  -1, 
b  *  -5,  c21  “  10,  and  the  other  c^j  are  zero. 

We  compare  the  boundary  conditions 


(3.17) 


v  -  0 
u  ■  g 


(x  -  0) 


and 

3  ac 

(3.18)  jj-j*  +  (g-~ ■")  u  -  0  (x  -  0)  . 

u  -  g 

Here  g  is  a  given  function  of  t.  The  simpler  conditions  (3.17)  are 
derived  from  the  method  of  characteristics  discussed  in  Section  3.1.  The 
first  condition  in  (3.18)  was  obtained  from  the  results  of  the  uncoupling 
process. 

In  our  computation  the  system  is  approximated  by  the  leap  frog 

difference  scheme.  The  function  g  in  the  boundary  conditions  is  equal 

to  a  half  period  of  a  sine  wave  which  is  extended  by  zero.  A  forward 

difference  is  used  to  approximate  the  derivative  in  (3.18).  The  surfaces 

/~2  2 

pictured  in  Figures  3.2  and  3.3  are  graphs  of  ru  +  v  as  a  function 
of  x  and  t.  In  Figure  3.1  we  illustrate  the  configuration  of  these 
surface  plots. 
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Figure  3.1 


In  the  computations  we  set  the  solution  equal  to  zero  when  t  *  0. 

The  nonzero  part  of  the  solution  is  thus  due  entirely  to  the  nonzero 
boundary  data,  so  it  is  possible  to  study  the  influence  of  the  boundary 
data  by  examining  the  size  of  the  solution  in  various  parts  of  the 
(x,t)  plane.  The  solution  corresponding  to  the  simple  boundary  condition 
(3.17)  is  graphed  in  Figure  3.2,  and  the  solution  corresponding  to  the 
more  refined  condition  (3.18)  is  given  in  Figure  3.3.  It  is  clear  from 
the  figures  that  the  second  condition  is  much  more  effective  in 
suppressing  the  fast  part  of  the  solution. 


Figure  3.2.  Solution  corresponding  to 


Figure  3. 


boundary  condition  (3.17) 


VARIABLE  COEFFICIENTS;  ANALYSIS  OF  THE  BOUNDARY  CONDITIONS 


In  this  section  we  discuss  a  couple  of  theoretical  questions  which  are 
suggested  by  the  calculations  given  in  Section  3. 

4.1  Generalization  to  Systems  Having  Variable  Coefficients. 

Problems  with  variable  coefficients  can  be  treated  by  using 
pseudo-differential  operators  in  a  manner  which  is  analogous  to  the 
development  in  the  preceding  section.  This  is  done  for  the 
one-dimensional  problem  in  [6] ,  and  in  Part  II  we  use  these  operators  on 
an  example  in  two  space  dimensions.  The  purpose  of  the  present  sub¬ 
section  is  to  give  a  brief  description  of  the  ideas  behind  the 
generalization. 

A  pseudo-differential  operator  has  the  form 
(4.1)  Pu(t)  -  /  ei^tp(t,5)  &(£)d€ 

For  example,  if  the  symbol  p  is  a  polynomial  in  €  with  coefficients 
depending  at  most  on  t,  then  P  is  also  a  differential  operator.  In 
the  Appendix  in  Part  II  we  describe  a  class  of  permissible  symbols  and 
summarize  some  formal  properties  of  these  operators. 

During  the  calculations  in  Section  3  we  multiplied  the  Fourier 
transform  by  functions  of  the  dual  variable  €.  In  effect,  we  were 
applying  pseudo-differential  operators  whose  symbols  were  independent  of 
t.  In  this  case  there  is  no  need  to  carry  out  the  integration  which 
inverts  the  Fourier  transform,  and  the  introduction  of  a  name  for  these 
operators  contributes  nothing  to  the  discussion. 

When  the  coefficients  in  the  system  depend  on  t,  the  manipulations 
in  Section  3  are  obviously  not  valid.  However,  it  is  possible  to  carry 
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out  analogous  calculations  which  are  valid  in  this  case,  instead  of 
multiplying  the  Fourier  transform  by  functions  of  the  dual  variable,  one 
applies  operators  defined  in  terms  of  the  inverse  Fourier  transform,  as  in 
(4.1).  These  calculations  involve  the  use  of  expansions  which  are  valid 
asymptotically  as  |€|  +  ",  and,  roughly  speaking,  the  leading  order 
terms  in  these  esqpansions  are  the  results  which  one  would  obtain  by 
calculating  formally  with  Fourier  transforms  as  in  Section  3.  This  is  to 
be  expected,  since  variable  coefficients  appear  almost  constant  to  a  wave 
whose  frequency  is  sufficiently  high.  The  expressions  for  the  lower  order 
terms  in  the  expansions  involve  derivatives  of  the  coefficients. 
Pseudo-differential  operators  can  thus  be  regarded  as  a  formalism  for 
studying  the  high-frequency  asymptotic  behavior  of  the  solution,  and  in 
the  present  work  they  are  used  both  to  uncouple  the  system  and  to  derive 
boundary  conditions  from  the  results  of  the  uncoupling. 

4. 2  Estimate  of  the  Size  of  the  Fast  Part  of  the  Solution 

Section  3.4  provides  numerical  evidence  of  the  effectiveness  of  the 
boundary  conditions  which  were  derived  in  Sections  3.2  and  3.3.  Here  we 
outline  a  method  for  justifying  the  boundary  conditions  analytically. 
Details  are  given  in  [6] .  This  method  also  works  for  problems  in  several 
space  dimensions. 

When  a  system  with  constant  coefficients  is  uncoupled  as  in  Section 
3,  one  can  obtain  an  ordinary  differential  equation  in  x  for  the  Fourier 
transform  of  the  portion  of  the  solution  which  is  entering  the  region  at 
the  fast  speed.  There  is  an  equation  for  each  frequency.  The  solution 
can  be  estimated  in  a  manner  which  is  essentially  the  same  as  the  one  used 
to  obtain  energy  estimates  for  hyperbolic  equations.  The  result  is  a 


bound  for  the  incoming  fast  part  in  terms  of  its  value  at  the  boundary  and 
the  magnitude  of  the  forcing  caused  by  other  components  in  the  system. 

The  choice  of  boundary  conditions  and  the  fact  that  the  forcing  is  small 
imply  that  the  incoming  fast  part  of  the  solution  must  be  small. 

When  the  coefficients  in  the  system  are  not  constant,  it  is  not 
possible  to  examine  the  solution  one  frequency  at  a  time.  However,  one 
can  still  consider  the  restriction  of  the  solution  to  a  conical 
neighborhood  of  a  bicharacteristic  curve.  This  amounts  to  a  study  of  the 
propagation  along  a  ray  of  a  family  of  neighboring  high  frequencies. 

Energy  estimates  for  the  restriction  yield  results  which  are  analogous  to 
those  found  in  the  case  of  constant  coefficients.  The  construction  of  a 
suitable  restriction  operator  is  given  in  Nirenberg  [10,  p.44]. 

In  several  dimensions  the  restriction  process  enables  one  to  examine 
the  effects  of  the  boundary  conditions  at  various  angles  of  incidence  to 
the  boundary.  The  utility  of  this  is  implied  by  the  discussion  appearing 


in  the  next  section 


5.  SYSTEMS  IN  SEVERAL  SPACE  DIMENSIONS 

In  this  section  we  generalize  the  methods  of  Section  3  to  systems  in 
several  dimensions.  The  main  point  of  discussion  is  the  problem  of  uncoupling 
the  leading  order  part  of  the  system.  The  lower  order  terms  can  be  treated  in 
the  same  manner  as  before.  A  summary  of  the  uncoupling  process  is  given  in 
the  last  sub-section. 


5. 1  Assumptions 

We  consider  the  hyperbolic  system 

(5.1)  wfc  =  Awx  +  BWy  +  Cw 

for  x  >  0,  y  €  R.  Here  w(x,y,t)  €  H P,  and  A,B,  and  C  are  real 
n  x  n  matrices.  We  assume  that  A  is  nonsingular,  i.e. ,  that  the 
boundary  is  noncharacteristic,  and  without  loss  of  generality  we  also 
assume  that  A  is  diagonal.  In  order  to  simplify  the  notation  we  have 
chosen  a  system  in  two  space  dimensions.  However,  it  will  be  clear  that 
the  discussion  is  equally  valid  for  systems  in  higher  dimensions  where 
x  >  0  and  y  €  for  k  >  2. 

There  is  no  serious  loss  of  generality  in  assuming  that  the  spatial 
domain  is  a  half-space.  If  the  given  domain  does  not  have  this  form  but 
still  has  a  smooth  boundary,  then  it  is  possible  to  localize  the  problem 
with  a  partition  of  unity  and  then  map  each  boundary  portion  into  the 
boundary  of  a  half-space.  In  the  new  coordinates  the  problem  will  have 
the  form  given  above. 


The  system  (5.1)  has  been  assumed  to  be  hyperbolic.  In  this  paper 


this  means  that  for  every  real  C  and  to  ,  the  symbol 

(S.2)  CA  +  (OB 

has  real  eigenvalues  and  a  complete  set  of  real  eigenvectors.  There  will 
be  no  need  to  assume  that  A  and  B  are  symmetric  or  that  the  system  is 
strictly  hyperbolic. 

In  order  to  have  a  system  with  at  least  two  time  scales,  we  assume 
that  certain  eigenvalues  of  the  symbol  (5.2)  are  substantially  greater  in 
magnitude  than  the  others.  Zn  the  case  of  the  linearized  shallow  water 
equations  this  symbol  has  eigenvalues  -u»cr  and  -u»a  +  c|a|,  where 
u  •  (u^Uj)  is  the  velocity  of  the  flow  about  which  the  system  has  been 
linearized,  and  o  is  the  vector  of  dual  variables  (C,<o).  If  |u|<<c, 
then  this  system  has  two  time  scales.  There  is  a  similar  set  of 
eigenvalues  for  the  three-dimensional,  five-component  Euler  system  for  gas 
dynamics.  Zn  this  case  the  small  eigenvalue  has  multiplicity  three. 

5. 2  An  Overview  of  the  Problem 

We  wish  to  transform  the  system  (5.1)  to  an  approximate  diagonal 
form,  or  at  least  block  diagonal  form,  so  that  each  of  the  new  dependent 
variables  can  be  identified  as  a  slow,  incoming  fast,  or  outgoing  fast 
portion  of  the  solution.  We  will  then  attempt  to  set  the  incoming  fast 
components  equal  to  zero. 

The  immediate  goal  is  to  uncouple  the  leading  order  terms  in  the 
system  (5.1).  After  this  has  been  accomplished,  we  can  use  the  methods  of 
Section  3  to  reduce  the  coupling  caused  by  the  lower  order  terms.  Zt 
would  suffice  to  obtain  a  block  diagonal  form  for  the  system,  instead  of 
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diagonal  form,  since  there  is  no  need  to  separate  various  incoming  fast 
components  or  various  incoming  or  outgoing  slow  components.  This 
situation  could  occur  with  the  Euler  equations,  for  example. 

He  have  assumed  that  the  matrix  A  in  (5.1)  is  already  in  diagonal 
form.  This  involves  no  loss  of  generality,  since  if  A  is  not  in  that 
form  we  can  find  a  similarity  transformation  which  makes  it  diagonal  and 
then  adopt  a  suitable  change  of  dependent  variable.  Unfortunately,  it  is 
not  true  in  general  that  this  transformation  can  also  uncouple  the 
matrix  B.  It  is  therefore  necessary  to  do  something  extra  if  we  want  to 
uncouple  the  entire  principal  part  of  (5.1). 

In  the  case  of  constant  coefficients  it  may  be  tempting  to  use 
Fourier  transforms  in  x  and  y.  This  would  yield  the  equation 

ft(C,<d,t)  -  (iCA  +  iWB)w  +  Cft. 

The  leading  symbol  of  this  equation  can  be  diagonalized  easily  because  it 
is  a  scalar  multiple  of  the  symbol  (5.2)  discussed  earlier.  However,  the 
use  of  Fourier  transforms  in  x  requires  the  use  of  information  about  the 
solution  away  from  the  boundary  x  -  0,  and  this  is  not  appropriate  in  a 
discussion  of  boundary  conditions.  He  will  instead  apply  Fourier 
transforms  in  time  and  in  the  tangent  variable  y.  The  use  of  the 
transform  in  t  can  be  justified  in  the  manner  indicated  in  Section  3.2. 

The  system  (5.1)  can  be  written  in  the  form 

(5.3)  wx  “  “  A_1Bwy  -  A-10 

Let  w(x,w,S)  denote  the  Fourier  transform  of  w  with  respect  to  y 
and  t  for  fixed  x.  Equation  (5.3)  implies 
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(5.4)  v x<xfw,5)  “  (Ua"1  -  i(*»A_1B)w  -  a"1#. 

It  is  necessary  to  determine  the  values  of  u  and  5  for  which  the 
symbol 

(5.5)  5a"1  -  «A-1B 

can  be  diagonalized,  or  brought  to  block  diagonal  form,  and  we  must 
determine  whether  such  an  uncoupling  can  produce  a  transformed  system  in 
which  each  component  of  the  dependent  variable  can  be  identified  as  slow, 
incoming  fast,  or  outgoing  fast.  The  answers  to  these  questions  are  not 
immediately  obvious,  since  we  are  applying  transforms  in  the  nonstandard 
variables  y  and  t. 

5.3  Properties  of  the  Symbol  (5.5) 

In  this  sub-section  we  consider  the  questions  posed  in  the  preceding 
paragraph,  interpret  the  results  physically,  and  discuss  methods  for 
obtaining  explicit  formulas  for  transformations  which  bring  about 
approximate  diagonal  forms  for  (5.5). 

We  must  first  consider  the  eigenvalues  and  eigenvectors  of  (5.5). 
Suppose  that  t  is  a  real  eigenvalue  and  that  v  is  a  corresponding 
eigenvector.  This  means  that 

(5.6)  (5a"1  -  oja"1B)v  =  Cv. 

If  we  multiply  by  A  and  rearrange  the  terms,  the  result  is 

(5.7)  (?A  +  UB)v  =  5v. 

The  matrix  CA  +  wb  is  the  symbol  (5.2)  which  we  would  obtain  by  writing 
the  system  in  the  more  common  form  (5.1)  and  then  applying  Fourier  trans- 
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forms  in  the  usual  variables  x  and  y.  According  to  (5.6)  and  (5.7), 
this  symbol  imposes  the  same  relations  between  the  dual  variables 
C,  w,  and  $  as  the  symbol  (5.5),  and  it  is  possible  to  find  the 
eigenvectors  of  one  symbol  by  examining  the  eigenvectors  of  the  other. 

The  difference  between  the  two  situations  is  that  in  one  case  the  variable 
C  is  treated  as  a  function  of  u>  and  £,  and  in  the  other  case  5  is 
treated  as  a  function  of  ?  and  <■>.  This  correspondence  between  the  two 
symbols  will  be  exploited  in  studying  (5.5),  since  at  this  point  in  the 
discussion  we  know  a  great  deal  more  about  (5.2)  than  we  do  about  (5.5). 

In  order  to  have  a  system  with  multiple  time  scales  we  have  assumed 
that  certain  eigenvalues  €  of  (5.2)  are  considerably  larger  than  the 
others.  An  example  of  such  a  set  of  eigenvalues  is  graphed  in  Figure 
5.1(a).  In  this  example  there  are  two  relatively  large  eigenvalues  and 
one  smaller  eigenvalue  for  each  C  and  u>.  This  is  the  configuration  for 
the  shallow  water  equations,  and  it  is  similar  to  the  configuration  for 
the  Euler  equations  of  gas  dynamics.  In  the  latter  case  the  small 
eigenvalue  has  multiplicity  three.  Throughout  this  discussion  we  will 
assume  that  the  largest  eigenvalues  of  (5.2)  have  graphs  which  are  narrow 
cones,  although  not  necessarily  right  circular  cones.  These  eigenvalues 
must  come  in  pairs,  since  if  (C,u,£)  is  a  solution  of  (5.7)  then  so  is 
(-5 ,-w,-5 ) .  These  large  eigenvalues  correspond  to  rapidly  moving  waves. 
The  fact  that  the  graphs  are  not  necessarily  right  circular  cones  means 
that  the  speed  can  vary  somewhat  with  the  direction  of  propagation.  We 
will  denote  by  Q  the  double  cone  which  corresponds  to  the  largest 
eigenvalues  of  (5.2),  and  we  will  denote  by  T  the  projection  of  ft  onto 
the  (u>,£)  space.  These  are  labeled  in  Figure  5.1. 


(a)  Graph  of  the  relations  (5.6),  (5.7)  (b)  Projection  onto  the 

(b>,£)  space 


Figure  5.1 

We  now  use  Figure  5.1(a)  to  study  the  behavior  of  the  eigenvalues 
C  of  the  symbol  (5.5),  ?A  1  -  u>a  1B.  First  of  all,  it  is  apparent  that 
the  number  of  real  eigenvalues  must  vary  with  the  position  of 
(u>,€).  if  (<■>,€)  lies  in  T,  then  there  are  two  real  values  of  C 
which  are  associated  with  the  surface  0.  As  ((■>,£)  approaches  the 
boundary  of  T,  these  values  of  t  coalesce,  and  when  (<•>,£)  leaves  T 
the  eigenvalues  leave  the  real  axis  and  form  a  pair  of  complex 
conjugates.  The  eigenvalues  cannot  be  real,  since  for  any  real  ?  the 
point  ((,<»,€)  must  lie  on  one  of  the  surfaces  in  Figure  5.1(a).  They 
are  complex  conjugates  because  they  are  eigenvalues  of  a  real  matrix. 

It  is  safe  to  assume  that  for  (w,S)  in  a  neighborhood  of  T  there 


is  no  problem  in  solving  for  the  values  of  C  associated  with  surfaces 


different  from  8.  In  the  case  of  the  shallow  water  equations  the  other 
£  satisfies  the  equation  5  =  -u^C  -  u^w.  We  have  assumed  that  the 
matrix  A  in  (5.1)  is  nonsingular,  which  in  this  case  is  equivalent  to 
saying  u^  ?  0.  It  is  therefore  possible  to  solve  for  £  in  terms  of 
oi  and  5,  whether  or  not  (u>,£)  is  in  I*.  A  similar  situation  holds 
for  the  Euler  equations  of  gas  dynamics. 

We  now  characterize  the  behavior  of  (5.5)  when  (<•>,£)  lies  in  T. 
It  will  become  clear  a  little  later  that  this  is  the  only  portion  of  the 
(<■>,£)  space  in  which  we  are  really  interested. 

Proposition.  If  (w,5)  is  in  T,  then  the  symbol  (5.5), 

Sa  1  -  oja  1B,  has  real  eigenvalues  and  a  complete  set  of  real 
eigenvectors.  This  is  not  the  case  if  ((■>,£)  is  not  in  I*.  The  eigen¬ 
vectors  can  be  determined  from  those  of  the  symbol  (5.2),  £A  +  u>B. 


Proof.  Equation  (5.6)  and  (5.7)  show  that  the  eigenvectors  of  (5.2) 
are  also  eigenvectors  of  (5.5).  We  know  that  (5.2)  has  a  complete  set  of 
real  eigenvectors  corresponding  to  fixed  (£,<*»)  and  various  eigenvalues 
€.  We  want  to  show  the  same  thing  for  (5.5),  for  fixed  (u,€)  in  T, 
and  various  eigenvalues  £ • 

Suppose  that  (<»>,£)  in  T,  and  let  £_,...,£  denote  the 

l  m 

eigenvalues  of  (5.5).  For  each  £^  choose  a  basis  Bj  for  the 
eigenspace  of  £^A  +  wB  corresponding  to  eigenvalue  £.  We  are  allowing 
for  the  possibility  that  the  symbol  (5.2)  might  have  mutliple 
eigenvalues.  The  elements  of  B  j  are  also  eigenvectors  of  (5.5) 
corresponding  to  the  eigenvalue  5^.  We  claim  that  the  union  of  the  Bj 
is  a  complete  set  of  vectors.  There  are  clearly  enough  of  these 
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vectors*  The  fact  that  they  are  linearly  independent  follows  from  an 
argument  which  is  essentially  the  one  which  shows  that  eigenvectors 
corresponding  to  distinct  eigenvalues  are  linearly  independent*  This 
completes  the  proof.  | 

The  matters  discussed  in  this  section  can  be  given  a  physical 
interpretation.  Suppose  that  the  coefficients  in  (5.1)  are  constant,  and 
let  C  «  0.  This  gives  the  system 

(5.8)  wt  ”  Awx  +  BWy. 

If  we  insert  a  plane  wave  solution  v  exp(i(x  +  iuty  +  i£t)  into  (5.8), 
where  v  is  a  ^ctor,  the  result  is  ?v  *  ((A  +  u>B)v.  This  is  the 
condition  (5.7)  which  was  discussed  earlier.  The  surfaces  in  graphs  like 
Figure  5.1(a)  therefore  define  the  set  of  all  possible  frequencies  for 
plane  wave  solutions  to  (5.8).  It  is  apparent  that  the  rapidly  moving 
waves  are  associated  with  ft,  which  is  why  we  are  interested  in  the 
behavior  of  (5.5)  only  for  (<*>,€)  in  I*.  This  analysis  is  also  valid 
when  C  +  0,  provided  we  are  considering  high-frequency  waves. 

In  graphs  like  Figure  5.1(a)  there  is  a  particular  wave  speed 
associated  with  each  surface  which  defines  C  as  a  function  of 
o>  and  £.  This  implies  that  it  is  possible  to  separate  fast  waves  from 
slow  waves  by  diagonalizing  the  symbol  (5.5).  It  is  also  necessary  to 
detect  the  directions  in  which  the  fast  waves  are  moving,  since  we  want 
the  diagonal  form  to  distinguish  incoming  and  outgoing  fast  waves.  This 
separation  is  possible.  The  vector  group  velocity  for  plane  wave 
solutions  of  (5.8)  is  equal  to  the  gradient  of  -5  with  respect  to 
(C,<*>).  From  this  it  follows  that  each  determination  of  ?  on  each 
section  of  T  can  be  identified  with  motion  into  or  out  of  the  domain 


x  >  0,  and  it  is  also  clear  that  the  edge  of  T  corresponds  to 

tangentially  moving  waves.  By  properly  defining  the  branches  of  C  on 

the  two  sections  of  T,  we  can  therefore  separate  the  fast  part  of  the 

solution  into  incoming  and  outgoing  components.  This  justifies  our 

decision  to  seek  diagonal  form  for  the  symbol  (5.5). 

He  need  to  say  a  little  more  about  the  directions  in  which  the 

various  waves  propagate.  When  we  seek  explicit  formulas  for  uncoupling 

the  symbol  (5.5)/  we  will  introduce  approximations  which  are  valid 

at  U) 

asymptotically  as  ^  ♦  0.  The  case  ^  *  0  corresponds  to  normal  phase 
velocity  for  a  wave  of  the  form  exp(i£x  +  iuy  +  i£t),  and  it  corresponds 
approximately  to  a  normal  group  velocity  for  fast  waves.  If  0  were  a 
right  circular  cone,  then  these  velocities  would  coincide.  The  approxi¬ 
mations  will  thus  lead  to  boundary  conditions  which  work  well  for  fast 
waves  traveling  in  directions  having  sizeable  normal  components,  but  they 
may  not  work  well  for  waves  moving  in  directions  which  are  nearly 
tangential,  i.e.,  for  («,£)  near  the  edge  of  T.  These  tangential  waves 
do  not  present  any  real  problem,  since  they  cannot  influence  the  interior 
very  rapidly.  The  approximation  schemes  are  therefore  worth  using. 

One  such  scheme  is  the  perturbation  method  given  in  the  Appendix. 

The  symbol  (5.5)  can  be  written  in  the  form  ?(A  1  -  ^A  ^B) •  When  |^| 

-1  <u  -1 

is  small,  the  matrix  A  -  ^A  B  is  a  perturbation  of  the  diagonal 
matrix  A-1,  and  the  perturbation  method  can  be  used  to  uncouple  it  to 

b> 

higher  powers  in  From  the  discussion  in  the  Appendix  it  is  apparent 

that  in  the  case  of  multiple  eigenvalues  this  method  cannot  guarantee 
diagonal  form,  but  instead  can  give  a  suitable  block  diagonal  form. 

Another  way  to  uncouple  (5.5)  is  to  compute  its  eigenvectors 
directly.  One  could  either  work  directly  with  (5.5)  or  instead  find 


i 


s 


i 


eigenvectors  of  the  symbol  (5.2),  ?A  +  wb,  and  then  use  the  ideas  of  the 

Proposition  to  translate  these  vectors  into  eigenvectors  of  (5.5).  This 

would  give  an  exact  diagonalization  of  (5.5)  when  (u,£)  is  in  T,  but 

the  egressions  for  the  eigenvectors  cam  be  complicated,  and  in  order  to 

obtain  local  boundary  conditions  it  would  usually  be  necessary  to 

approximate  these  vectors  with  simple  expressions.  We  would  again  use 

a) 

approximations  which  are  valid  asymptotically  as  ^  ♦  0.  Unlike  the 
perturbation  method,  which  employs  one  fixed  method  of  approximation,  this 
approach  allows  the  use  of  various  Taylor  or  rational  Padfe  approxi¬ 
mations.  Engquist  and  Majda  ([1],[2])  found  the  latter  approximations 
particularly  useful  in  their  work  on  absorbing  boundary  conditions  for 
scalar  wave  equations.  The  greater  flexibility  of  this  approach  may  be  an 
advantage,  in  general,  but  in  the  example  appearing  in  Part  II  we  are  able 
to  use  the  perturbation  method  and  obtain  satisfactory  results. 

5.4  Summary  of  the  Uncoupling  Process 

We  now  outline  the  uncoupling  process  in  terms  of  formal 
manipulations  of  Fourier  transforms.  In  Part  II  the  method  will  be 
translated  into  the  language  of  pseudo-differential  operators  and  applied 
in  detail  to  the  linearized  shallow  water  equations. 

We  first  solve  for  wx  in  (5.1)  to  obtain 
wx  *  A~^wt  -  A-1Bwy  -  A^Cw.  In  order  to  simplify  the  notation  we  change 
the  meaning  of  A,  B,  and  C  and  write  the  system  as 

(5.9)  wx  "  Awt  +  BWy  +  Cw 

The  matrix  A  is  diagonal. 
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When  w  apply  Fourier  transforms  in  y  and  t  to  (5.9),  the  result  Is 

(5.10)  ft  -  (15a  +  iu>B)w  +  Cw 

The  symbol  £ A  +  wb  Is  the  same  as  the  symbol  (5.5)  discussed  earlier, 
aside  from  the  change  in  notation.  According  to  the  remarks  of  the 
preceding  sub-section,  this  symbol  is  diagonalizable  for  (<*>,£)  in  T, 
and  it  is  only  for  («,£)  in  T  that  we  can  have  rapidly  moving  waves. 
For  the  sake  of  neatness  we  use  a  cutoff  function  to  restrict  attention  to 
that  set.  Let  +  be  a  smooth  function  of  w  and  ?  which  is  equal  to 
zero  outside  T  and  which  is  equal  to  1  on  all  of  T  except  for  a  thin 
layer  near  the  boundary.  Equation  (5.10)  can  then  be  written 

(5.11)  wx(x,w,€)  -  (i?A  +  iu$B)w  +  Cw  +  ic«»(  1  -  +)Bw 

The  last  term  in  (5.11)  is  an  error  term  which  we  will  denote  by  Ew.  It 
is  nonzero  only  near  the  edge  of  T,  and  for  fast  waves  this  corresponds 
to  nearly  tangential  incidence.  The  error  term  is  therefore 
insignificant. 

The  next  task  is  to  uncouple  the  leading  order  terms.  That  is,  we 
find  a  matrix  q  such  that 

(5.12)  q(i€A  +  iw$B)q-1 

(A) 

is  closer  to  diagonal  form,  or  to  a  suitable  block  diagonal  form,  when 

is  small.  This  question  was  discussed  in  Section  5.3.  Denote  (5.12)  by 

g  +  r,  where  g  is  uncoupled  and  r  is  the  error  in  the  uncoupling. 

3  “1 

Equation  (5.11)  then  becomes  r— (qw)  ■  q(  i?A  +  iu>$B)q  qw  + 

0  x 


(5.13) 


1 

I 

1 

B 


3C, 


-  gwo  +  (qC  +  qx)q'1wo  + 


o  a. 


where  w  ■  qw  and  E  ■  r  +  qEq  ^  *  E  w  is  an  error  term 
o  o  o  o 

(l) 

which  ie  snail  when  is  close  to  zero. 

Before  we  proceed  further,  it  is  necessary  to  discuss  the  degree  of 

homogeneity  in  »  and  £,  or  "order",  of  various  expressions  which 

appear  in  (5.13).  The  matrix  q  can  be  assumed  to  be  homogeneous  of 

degree  zero  for  large  u  and  £,  since  we  can  first  uncouple  the  symbol 

2  2 

i£A  +  iw$B  when  #  +  (  ■  1  and  then  extend  q  to  be  constant  along 

rays  from  the  origin.  In  the  example  in  Part  II  we  will  use  one  iteration 
of  the  perturbation  method,  so  that  q  will  have  the  form  q  *  I  +  ^M. 

(i)  2 

In  this  case  the  error  in  the  uncoupling  will  be  r  ■  0(€(g4)  ). 

Because  q  has  order  zero,  (5.12)  implies  that  the  uncoupled  symbol  g 
must  have  order  one.  Similarly,  the  product  (qC  +  qx)q  '  must  have 
order  zero.  The  system  (5.13)  is  therefore  uncoupled  to  order  zero. 


modulo  the  insignificant  error  term  Eqwq. 


For  simplicity  we  will  let  z  ”  (qC  +qx)q-1  and  write  (5.13)  as 

3  w 

(5.14) 


3w 
_ o 

3x 


gw  +  zw  +  E  w 
o  o  o  o 


When  the  coefficients  in  the  system  depend  on  y  or  t,  the  coupling  of 
order  zero  must  also  include  some  extra  terms  generated  by  the  uncoupling 
of  the  leading  order  part  of  the  system.  These  extra  terms  will  appear  in 
the  discussion  in  Part  II. 

We  now  use  the  method  of  Section  3. 2  to  uncouple  the  terms  of  order 
zero.  We  first  multiply  (5.14)  by  a  matrix  of  the  form  I  +  k,  where 
k  is  a  matrix  of  order  -1  which  is  to  be  determined,  to  obtain 
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(5.15) 


t|[(I  +  k)wl  -  (I  +  k)g(I  +  k)”1 (I  +  k)* 

0X0  o 

+  (I  +k)zw  +  (I  +  k)E  w  +  k  w  . 

o  oo  xo 

Let  w  -  (I  +  k)w  ,  E  -  (I  +  k)E  (I  +  k)”1,  and  write 
1  o  1  o 

(I  +  k)-1  as  I  -  k  +  k2  ..  .  The  system  (5.15)  becomes 
3  w 

—3-^  ■  (I  +  k)g(I  -  k)w  +  zw  +  order(-1)w  +  E  w 

OX  II  111 

*  gw^  +  (kg  -  gk  +  z)w^  +  orderf-Dw^  +  E^ 

To  eliminate  the  coupling  o£  order  zero,  we  choose  k  so  that 

kg  -  gh  +  z  is  diagonal,  or  block  diagonal.  The  matrix  k  will  in  fact 

have  order  -1,  as  we  assumed  above,  since  g  has  order  one  and  z  has 

order  zero.  As  noted  in  Section  3.2,  this  method  can  be  repeated  to 

reduce  further  the  coupling  caused  by  lower  order  terms. 

The  new  dependent  variable  is  defined  by 
-  (I+k)wo»  (1+  k)qw.  To  obtain  local  boundary  conditions,  we 
express  the  incoming  fast  components  of  w^  in  terms  of  w,  set  these 
components  equal  to  zero  at  the  boundary,  clear  denominators,  and  invert 
the  Fourier  transform.  An  analogue  for  systems  with  variable  coefficients 


will  be  given  in  Part  II 


APPENDIX.  A  Perturbation  Lemma 

Here  we  present  a  method  for  reducing  the  coupling  found  in  matrices  which 
are  perturbations  of  block  diagonal  matrices.  This  method  can  be  used  to 
partially  uncouple  the  leading  symbol  in  the  system  (5.11),  and  it  is 
essentially  the  method  which  has  been  used  to  reduce  the  coupling  caused  by 
lower  order  terms.  We  present  it  as  a  separate  lemma  for  the  sake  of  clarity 
and  generality.  Various  versions  of  this  method  have  been  used  in  [8],  [11], 
and  [12]. 

Proposition.  Let  A  and  B  be  square  matrices  of  equal  dimension. 
Suppose  that  A  is  block  diagonal,  and  let  A,,,...,An  denote  the  blocks 
on  the  diagonal.  If  no  two  of  the  Aj  have  any  eigenvalues  in  common, 
then  for  small  €  the  sum  A  +  €B  can  be  uncoupled  to  order  e2.  More 
precisely,  there  exists  a  matrix  M  such  that  for  €  sufficiently  small, 

(I  +  ®1)(A  +  eBHI  +  8M)"1  -  A  +  e*  (block  diagonal  matrix)  +  0(e2). 

A  method  for  constructing  M  will  be  given  in  the  proof. 

Proof .  For  small  e  the  inverse  (I  +  €M)_1  exists  and  is  equal  to 
I  -  04  +  e^2  -...  .  We  can  therefore  write 

(I  +  0*)(A  +  €B) (I  +  €M)-1  »  (I  +  €M) (A  +  €B)(I  -  +  0(62)) 

*  A  +  €(MA  -  AM  +B)  +  0(€2). 

Our  goal  is  to  choose  M  so  that 

(1)  MA  -  AM  +  B 


is  block  diagonal.  For  the  sake  of  notation  we  will  partition  M  and  B 


into  block  structures  which  match  the  block  structure  of  A.  M 


and 


j  will  denote  the  blocks  in  the  (i,j)  position.  They  are  not 
necessarily  square,  since  we  are  not  assuming  that  A^  and  Aj  have  the 
same  dimensions.  The  (i,j)  block  in  (1)  can  then  be  written  as 
M^Aj  -  A^M^  +  j .  For  i  ^  j  we  want  this  to  be  equal  to  zero. 

He  are  therefore  faced  with  the  problem  of  solving  the  equation 


\ 

*  t 

i 

! 

1 


(2)  ^i jAj  "  AiMij  *  "Bij 

for  Once  we  have  done  this,  the  proof  is  complete.  There  are  no 

conditions  imposed  on  the  diagonal  blocks  M^,  so  these  may  be  chosen 
arbitrarily. 

If  A^  and  Aj  are  both  1x1  matrices,  i.e.,  scalars,  then  we 
obviously  need  to  have  A^  /  A j  in  order  to  be  able  to  solve  ( 2)  for 
arbitrary  B^j.  In  the  general  case  the  system  (2)  is  solvable  if  and 
only  if  A^  and  Aj  have  disjoint  spectra.  Proofs  of  this  fact  can  be 
found  in  several  different  references.  We  give  one  here  for  the  sake  of 
completeness • 

In  order  to  simplify  the  notation  we  will  write  (2)  in  the  form 

XS  -  TX  ■  Y,  where  S,  T,  and  Y  are  given  and  X  is  to  be 

determined.  We  are  assuming  that  S  and  T  are  square  matrices  which  do 
not  have  any  eigenvalues  in  common.  There  is  no  need  to  assume  that  they 
have  the  same  dimension.  We  will  denote  the  columns  of  X  and  Y  by  Xj 

and  yj,  and  we  will  denote  the  entries  of  S  by  s^j. 

We  can  assume  that  S  is  upper  triangular,  since  otherwise  we  can 
use  a  similarity  transformation  to  reduce  the  problem  to  that  case.  We 
will  solve  for  the  columns  of  X,  starting  from  the  left.  We  first 


1 
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have  s^x.,  -  Tx1  -  y^  The  matrix  s^Z  -  T  is  nonsingular  since 
Sjj  is  an  eigenvalue  of  S  and  therefore  not  an  eigenvalue  of  T.  The 
column  x^  is  therefore  determined  uniquely.  We  next  have 
s22x2  ~  Tx2  ~  ~  B^2X^’  T^is  system  has  a  unique  solution 

x2  since  s22  is  not  an  eigenvalue  of  T.  We  can  continue  in  this 
manner  to  solve  for  all  of  X.  We  note  that  the  condition  on  the 
eigenvalues  of  S  and  T  is  necessary  as  well  as  sufficient ,  since  it  is 
equivalent  to  the  statement  that  s^I  -  T  is  nonsingular  for  all  i. 

This  completes  the  proof  of  the  lemma ,  and  therefore  the  proof  of  the  main 
proposition. 


This  method  can  be  applied  repeatedly  to  uncouple  the  given  matrix  to 

2 

any  order.  To  eliminate  the  coupling  of  order  e  ,  we  would  introduce  a 
matrix  of  the  form  I  +  €^M2  an<*  solve  for  M2  in  the  manner  indicated 
above,  in  general,  to  reduce  the  coupling  from  order  e11  to  order 
€°+1,  we  would  use  a  transformation  of  the  form  I  + 
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